Event History Analysis II — Cox proportional hazard regression

Dagen program

  • Øvelsen til forelæsning 3.

  • Data og data management

  • Cox-regression

  • Analyse i praksis

  • Øvelse

Forventet udbytte

Viden

  • Cox-regressionen og dens forudsætninger.
  • Cox-regressionens fordele (og ulemper).

Færdigheder

  • Kodning af tidsvarierende variable til tidskonstante variable.
  • Grafisk præsentation af resultater fra en cox-regression.

Kompetancer

  • Analyse af Cox-regression i R.

Subject-Period file og single-episode file

Hvor vi sidste gang arbejde med en single-episode file, arbejder vi i konteksten af regressionsanalyse med mikrodata, ofte med en Subject-Period file — dvs. en fil, hvor hvert individ fremgår flere gange over tid.

Konstruktion af event file

Den centrale kode-mæssige udfordring er at konstruere event variablen.

03:00

klik for at se koden 
library(tidyverse)
library(readxl)
library(gt)

read_excel("eksempel_data.xlsx") %>% 
  gt(groupname_col = "grp") %>% 
  tab_header(
    title = "'Rå' data",
    subtitle = glue::glue("Fiktive data")
  ) %>% 
  tab_options(table.width = pct(40),
              table.font.size = figur_tekst_str) %>% 
  opt_table_lines("none")
'Rå' data
Fiktive data
id år børn alder region18 tid event
1
1 2000 0 18 Nord 1 0
1 2001 0 19 Nord 2 0
1 2002 1 20 Syd 3 1
1 2003 1 21 Syd 4 0
1 2004 1 22 Syd 5 0
2
2 2003 0 20 Syd 3 0
2 2004 0 21 Øst 4 0
2 2005 0 22 Øst 5 0
2 2006 0 23 Nord 6 0
2 2007 0 24 Nord 7 0
3
3 1999 0 18 Nord 1 0
3 2000 0 19 Nord 2 0
3 2001 0 20 Vest 3 0
3 2002 0 21 Vest 4 0
3 2003 1 22 Nord 5 1
4
4 1993 0 23 Øst 6 0
4 1994 0 24 Øst 7 0
4 1995 0 25 Øst 8 0
4 1996 1 26 Øst 9 1
4 1997 1 27 Øst 10 0

Konstruktion af event file

klik for at se koden 
library(tidyverse)
library(readxl)
library(gt)

read_excel("eksempel_data.xlsx") %>% 
  gt(groupname_col = "grp") %>% 
  tab_header(
    title = "'Rå' data",
    subtitle = glue::glue("Fiktive data")
  ) %>% 
  tab_options(table.width = pct(40),
              table.font.size = figur_tekst_str) %>% 
  tab_style(
    style = list(
      cell_fill(color = "pink")
      ),
    locations = cells_body(
      columns = event,
      rows = event == 0 & børn == 1
    )
  ) %>% 
  tab_style(
    style = list(
      cell_fill(color = "wheat")
      ),
    locations = cells_body(
      columns = event,
      rows = event == 1 & børn == 1
    )
  ) %>% 
  tab_style(
    style = list(
      cell_fill(color = "steelblue")
      ),
    locations = cells_body(
      columns = event,
      rows = event == 0 & år == 2007
    )
  ) %>% 
  opt_table_lines("none")
'Rå' data
Fiktive data
id år børn alder region18 tid event
1
1 2000 0 18 Nord 1 0
1 2001 0 19 Nord 2 0
1 2002 1 20 Syd 3 1
1 2003 1 21 Syd 4 0
1 2004 1 22 Syd 5 0
2
2 2003 0 20 Syd 3 0
2 2004 0 21 Øst 4 0
2 2005 0 22 Øst 5 0
2 2006 0 23 Nord 6 0
2 2007 0 24 Nord 7 0
3
3 1999 0 18 Nord 1 0
3 2000 0 19 Nord 2 0
3 2001 0 20 Vest 3 0
3 2002 0 21 Vest 4 0
3 2003 1 22 Nord 5 1
4
4 1993 0 23 Øst 6 0
4 1994 0 24 Øst 7 0
4 1995 0 25 Øst 8 0
4 1996 1 26 Øst 9 1
4 1997 1 27 Øst 10 0

Konstruktion af event file: Et eksempel

Funktionerne lag() og lead() er centrale, når vi skal konstruere vores event variabel. Disse funktioner bruges til at hente data fra henholdsvis den foregående og den efterfølgende (observations)række i vores dataframe. Altså, lag() og lead() bruges til at “kalde” oplysninger fra henholdsvis før og efter en given observationsrække.

Spg.: Hvordan løser det vores “problem” med at skabe vores event indikator?

02:00

Indikator for startstidspunktet per person:

library(tidyverse)

df.rå <- read_excel("eksempel_data.xlsx")

df.lagged <- 
  df.rå %>% 
  mutate(person_start = case_when(lag(id) == id ~ 0, TRUE ~ 1))
Konstruér starttidspunkt i forløbet
Fiktive data
id år børn alder region18 tid event person_start
1
1 2000 0 18 Nord 1 0 1
1 2001 0 19 Nord 2 0 0
1 2002 1 20 Syd 3 1 0
1 2003 1 21 Syd 4 0 0
1 2004 1 22 Syd 5 0 0
2
2 2003 0 20 Syd 3 0 1
2 2004 0 21 Øst 4 0 0
2 2005 0 22 Øst 5 0 0
2 2006 0 23 Nord 6 0 0
2 2007 0 24 Nord 7 0 0
3
3 1999 0 18 Nord 1 0 1
3 2000 0 19 Nord 2 0 0
3 2001 0 20 Vest 3 0 0
3 2002 0 21 Vest 4 0 0
3 2003 1 22 Nord 5 1 0
4
4 1993 0 23 Øst 6 0 1
4 1994 0 24 Øst 7 0 0
4 1995 0 25 Øst 8 0 0
4 1996 1 26 Øst 9 1 0
4 1997 1 27 Øst 10 0 0

Indikator for sluttidspunktet per person:

library(tidyverse)

df.rå <- read_excel("eksempel_data.xlsx")

df.lagged <- 
  df.rå %>% 
  mutate(person_start = case_when(lag(id) == id ~ 0, TRUE ~ 1)) %>% 
  mutate(person_slut = case_when(lead(id) == id ~ 0, TRUE ~ 1))
Konstruér starttidspunkt i forløbet
Fiktive data
id år børn alder region18 tid event person_start person_slut
1
1 2000 0 18 Nord 1 0 1 0
1 2001 0 19 Nord 2 0 0 0
1 2002 1 20 Syd 3 1 0 0
1 2003 1 21 Syd 4 0 0 0
1 2004 1 22 Syd 5 0 0 1
2
2 2003 0 20 Syd 3 0 1 0
2 2004 0 21 Øst 4 0 0 0
2 2005 0 22 Øst 5 0 0 0
2 2006 0 23 Nord 6 0 0 0
2 2007 0 24 Nord 7 0 0 1
3
3 1999 0 18 Nord 1 0 1 0
3 2000 0 19 Nord 2 0 0 0
3 2001 0 20 Vest 3 0 0 0
3 2002 0 21 Vest 4 0 0 0
3 2003 1 22 Nord 5 1 0 1
4
4 1993 0 23 Øst 6 0 1 0
4 1994 0 24 Øst 7 0 0 0
4 1995 0 25 Øst 8 0 0 0
4 1996 1 26 Øst 9 1 0 0
4 1997 1 27 Øst 10 0 0 1

Hvis vi nu vil definere begivenheden “fødsel af første barn” er udfordringen at “lokalisere” der hvor kvinden går fra 0 børn til 1 barn.

Jeg bruger her både case_when() og if_else() for illustrationen til begge funktinoer.

library(tidyverse)

df.rå <- read_excel("eksempel_data.xlsx")

df.lagged <- 
  df.rå %>% 
  mutate(person_start = case_when(lag(id) == id ~ 0, TRUE ~ 1)) %>% 
  mutate(person_slut = case_when(lead(id) == id ~ 0, TRUE ~ 1)) %>% 
  group_by() %>% 
  mutate(begivenhed = if_else(lag(børn) == 0 & børn == 1 & person_start == 0, 1, 0)) %>% 
  ungroup()
Konstruér indikator for begivenhed
Fiktive data
id år børn alder region18 tid event person_start person_slut begivenhed
1
1 2000 0 18 Nord 1 0 1 0 0
1 2001 0 19 Nord 2 0 0 0 0
1 2002 1 20 Syd 3 1 0 0 1
1 2003 1 21 Syd 4 0 0 0 0
1 2004 1 22 Syd 5 0 0 1 0
2
2 2003 0 20 Syd 3 0 1 0 0
2 2004 0 21 Øst 4 0 0 0 0
2 2005 0 22 Øst 5 0 0 0 0
2 2006 0 23 Nord 6 0 0 0 0
2 2007 0 24 Nord 7 0 0 1 0
3
3 1999 0 18 Nord 1 0 1 0 0
3 2000 0 19 Nord 2 0 0 0 0
3 2001 0 20 Vest 3 0 0 0 0
3 2002 0 21 Vest 4 0 0 0 0
3 2003 1 22 Nord 5 1 0 1 1
4
4 1993 0 23 Øst 6 0 1 0 0
4 1994 0 24 Øst 7 0 0 0 0
4 1995 0 25 Øst 8 0 0 0 0
4 1996 1 26 Øst 9 1 0 0 1
4 1997 1 27 Øst 10 0 0 1 0 
library(tidyverse)

df.rå <- read_excel("eksempel_data.xlsx")

df.lagged <- 
  df.rå %>% 
  mutate(person_start = case_when(lag(id) == id ~ 0, TRUE ~ 1),
         person_slut = case_when(lead(id) == id ~ 0, TRUE ~ 1)) %>% 
  group_by(id) %>% 
  mutate(begivenhed = if_else(lag(børn) == 0 & børn == 1 & person_start == 0, 1, 0),
         periode_slut = case_when(begivenhed == 1 | børn == 0 & person_slut == 1 ~ 1, TRUE ~ 0)) %>% 
  slice(1:which.max(periode_slut == 1)) %>% 
  ungroup()
Konstruér subject-period file
Fiktive data
id år børn alder region18 tid event person_start person_slut begivenhed periode_slut
1
1 2000 0 18 Nord 1 0 1 0 0 0
1 2001 0 19 Nord 2 0 0 0 0 0
1 2002 1 20 Syd 3 1 0 0 1 1
2
2 2003 0 20 Syd 3 0 1 0 0 0
2 2004 0 21 Øst 4 0 0 0 0 0
2 2005 0 22 Øst 5 0 0 0 0 0
2 2006 0 23 Nord 6 0 0 0 0 0
2 2007 0 24 Nord 7 0 0 1 0 1
3
3 1999 0 18 Nord 1 0 1 0 0 0
3 2000 0 19 Nord 2 0 0 0 0 0
3 2001 0 20 Vest 3 0 0 0 0 0
3 2002 0 21 Vest 4 0 0 0 0 0
3 2003 1 22 Nord 5 1 0 1 1 1
4
4 1993 0 23 Øst 6 0 1 0 0 0
4 1994 0 24 Øst 7 0 0 0 0 0
4 1995 0 25 Øst 8 0 0 0 0 0
4 1996 1 26 Øst 9 1 0 0 1 1

Herfra er en Single-episode file ligefrem at konstruere, ved:

library(tidyverse)

df.rå <- read_excel("eksempel_data.xlsx")

df.lagged <- 
  df.rå %>% 
  mutate(person_start = case_when(lag(id) == id ~ 0, TRUE ~ 1)) %>% 
  mutate(person_slut = case_when(lead(id) == id ~ 0, TRUE ~ 1)) %>% 
  group_by(id) %>% 
  mutate(begivenhed = if_else(lag(børn) == 0 & børn == 1 & person_start == 0, 1, 0)) %>% 
  mutate(periode_slut = case_when(begivenhed == 1 | børn == 0 & person_slut == 1 ~ 1, TRUE ~ 0)) %>% 
  ungroup() %>% 
  filter(periode_slut == 1)
Konstruér single-episode file
Fiktive data
grp id år børn alder region18 tid event person_start person_slut begivenhed periode_slut
1 1 2002 1 20 Syd 3 1 0 0 1 1
2 2 2007 0 24 Nord 7 0 0 1 0 1
3 3 2003 1 22 Nord 5 1 0 1 1 1
4 4 1996 1 26 Øst 9 1 0 0 1 1

Konstruktion af event file

Hvad så med en Single-episode file i stil med First_child.rda. Hvad er fordele og begrænsninger ved denne?

05:00

klik for at se koden 
df.rå <- read_excel("eksempel_data.xlsx")


df.lagged <- 
  df.rå %>% 
  mutate(person_start = case_when(lag(id) == id ~ 0, TRUE ~ 1)) %>% 
  mutate(person_slut = case_when(lead(id) == id ~ 0, TRUE ~ 1)) %>% 
  group_by(id) %>% 
  mutate(begivenhed = if_else(lag(børn) == 0 & børn == 1 & person_start == 0, 1, 0)) %>% 
  mutate(periode_slut = case_when(begivenhed == 1 | børn == 0 & person_slut == 1 ~ 1, TRUE ~ 0)) %>% 
  ungroup() %>% 
  filter(periode_slut == 1)


df.lagged %>% 
  gt() %>% 
  tab_header(
    title = "Konstruér single-episode file",
    subtitle = "Fiktive data"
  ) %>% 
  tab_options(table.width = pct(40),
              table.font.size = figur_tekst_str) %>% 
  opt_table_lines("none") %>% 
  tab_style(
    style = list(
      cell_fill(color = "wheat")
      ),
    locations = cells_body(
      columns = c(person_start, person_slut, periode_slut)
    )
  ) %>% 
  tab_style(
    style = list(
      cell_fill(color = "steelblue")
      ),
    locations = cells_body(
      columns = begivenhed
    )
  ) %>% 
  tab_style(
    style = list(
      cell_fill(color = "lightblue")
      ),
    locations = cells_body(
      columns = begivenhed,
      rows = event == 1
    )
  ) %>% 
  tab_style(
    style = list(
      cell_fill(color = "pink")
      ),
    locations = cells_body(
      columns = event,
      rows = event == 1
    )
  )
Konstruér single-episode file
Fiktive data
grp id år børn alder region18 tid event person_start person_slut begivenhed periode_slut
1 1 2002 1 20 Syd 3 1 0 0 1 1
2 2 2007 0 24 Nord 7 0 0 1 0 1
3 3 2003 1 22 Nord 5 1 0 1 1 1
4 4 1996 1 26 Øst 9 1 0 0 1 1

Cox proportional hazard model

Cox PH: relevans og popularitet

  • Cox PH-modellen er semi-parametrisk, hvilket betyder, at vi ikke behøver at specificere den underliggende sandsynlighedstæthedsfunktion (pdf) for hændelsestiden. (Genbesøg noter og slides til lektion 2, for introduktion til relevante pdf.) Dette gør modellen fleksibel, da vi kun estimerer effekten af kovariaterne på hazardraten, mens baseline hazard forbliver ukendt.

  • Modellen er robust og kan anvendes til data. Uanset om den underliggende fordeling for hændelsestiderne er Weibull, eksponentiel, lognormal eller en anden form, vil Cox-modellen give gode approksimationer, så længe den proportionale hazard-antagelse er opfyldt.

    • En væsentlig forudsætning for modellen er, at hazardraten for de forskellige grupper er proportionale over tid. Denne antagelse skal testes, da overtrædelser kan påvirke modellens anvendelighed.

  • En af de primøre fordele ved Cox-modellen er, at den nemt kan udvides til at inkludere tidsvarierende kovariater.

Cox proportional hazard model med fixed covariater.

For den \(i\)’te person til tiden \(t\) specificeres hazardfunktionen som:

\[ h_{i}(t)=h_{0}(t) \times \exp(\beta_{1} x_{1} + \dots + \beta_{k} x_{k}) \]

hvor,

  • \(h_{0}(t)\) (baseline hazard) angiver den underliggende risiko ved tid tt for en referenceperson, dvs. for en person med \(x_{1,i} = x_{2,i} = \dots = x_{k,i} = 0\). \(h_{0}(t)\) specificeres ikke a priori, hvilket betyder, at vi ikke antager nogen konkret form for dens fordeling.

  • \(\exp(\beta_{1} x_{1} + \dots + \beta_{k} x_{k})\) er den multiplicative effekt af de \(k\) covariater på hazardraten.

  • Bemærk at der ikke er noget egentlig intercept, da den faste effekt, som et intercept normalt ville repræsentere, er integreret i \(h_{0}(t)\). Ved estimering (via partial likelihood) udtrækkes koefficienterne \(\beta\), mens \(h_{0}(t)\) forbliver uspecificeret.

Cox proportional hazard model med fixed covariater.

For den \(i\)’te person til tiden \(t\) specificeres hazardfunktionen som:

\[ h_{i}(t)=h_{0}(t) \times \exp(\beta_{1} x_{1} + \dots + \beta_{k} x_{k}) \]

Dermed kan vi udtrykke hazardfunktionen for en covariatvektor \(\mathbf{X}\) som:

\[ h(t|\mathbf{X}) = h_{0}(t) \exp(\beta^{T}\mathbf{X}) \]

Cox proportional hazard model med fixed covariater. Sammenligning af grupper.

Med afsæt i data fra sidste øvelse, First_child.rda, antager vi, at vi har en gruppeindikator \(x\), som defineres som en dummyvariabel:

\[ x = \left\{\begin{matrix} 0, \: \text{hvis personen tilhører referencegruppen} \:\:\:\:\:\:\:\: \\ 1, \: \text{hvis personen ikke tilhører referencegruppen} \end{matrix}\right. \]

For personer i referencegruppen (\(x=0\)):

\[ h_{i}(t)=h_{0}(t) \times \exp(\beta \times 0) = h_{0}(t) \]

For personer i den sammenlignede gruppe (\(x=1\)):

\[ h_{i}(t)=h_{0}(t) \times \exp(\beta \times 1) = h_{0}(t) \times \exp(\beta) \]

Cox proportional hazard model med fixed covariater. Hazard-ratioen

Hazard-ratioen (HR) udtrykker forholdet mellem hazardfunktionen for en person uden for referencegruppen og en person i referencegruppen:

\[ \require{cancel} HR = \frac{h_{i}(t)}{h_{j}(t)} = \frac{ h_{0}(t) \times \text{exp}(\beta)}{h_{0}(t)} = \frac{ \cancel{h_{0}(t)} \times \text{exp}(\beta)}{\cancel{h_{0}(t)}} = \text{exp}(\beta) \]

Konceptuelt svarer HR til Odds Ratio (OR) i logistisk regression, idet den udtrykker den relative ændring i hazardraten mellem grupper:

  • Hvis HR er \(=1\) har personen uden for referencegruppen samme hazard-rate som personen i referencegruppen.

  • Hvis HR er \(<1\) har personen uden for referencegruppen en lavere hazard-rate end personen i referencegruppen.

  • Hvis HR er \(>1\) har ersonen uden for referencegruppen en højere hazard-rate end personen i referencegruppen.

Den proportionale hazard model

Fordi \(h_{0}(t)\) bliver ophævet i tæller og nævner, gælder det at hazardrater mellem to personer ved samme tid er uafhængig af selve tiden. Det skyldes, at baseline hazardfunktionen \(h_{0}(t)\) ophæves, når vi danner forholdet (hazard-ratioen) mellem to individer. Dermed er hazard-ratioen udelukkende bestemt af de eksponentierede kovariat-effekter og antages at være konstant over tid.

Den proportionale hazard model

Fordi hazard-ratioen (HR) mellem to personer er konstant over tid, gælder det at:

\[ \hat{HR} = \frac{h_{i}(t)}{h_{j}(t)} = \theta \] hvor \(\theta\) er en konstant, der udelukkende bestemmes af de eksponentierede kovariateffekter.

Den proportionale hazard model

Beskrevet andeledes: alle starter med den samme “basisrisiko” for at blive ramt af en begivenhed. Når vi sammenligner to personer, ser vi kun på, hvordan deres individuelle egenskaber gør, at deres risikoer afviger fra denne “basisrisiko” – og det forhold (hazard-ratioen) er konstant over tid.

Fx: Gruppe 1 har en højere risiko; og denne risiko er \(n\) højere end gruppe 2 over hele forløbet.

03:00

adjusted overlevelseskurver

Den centrale interesse er fortsat på at estimere overlevelseskurven, \(S(t)\). Som I ved fra de forrige lektioner er der en relation mellem \(S(t)\) og \(h(t)\) (genlæs noter eller slide til lektion 2+3).

Når vi inkluderer individuelle karakteristika (covariater) i analysen, udtrykker vi hazardfunktionen for en person med covariatvektoren \(\mathbf{X}\) som: \[ h(t, \mathbf{X}) = h_{j}(t) \times \text{exp}(\mathbf{X}\beta) \]

hvor

  • \(h_{j}(t)\) er baseline hazard: den risiko, der gælder for referencegruppen

  • \(\text{exp}(\mathbf{X}\beta)\) beskriver, hvordan de enkelte covariater ændrer denne risiko.

adjusted overlevelseskurver

Dette giver os en “kontrolleret” overlevelsesfunktion, som justeres for covariaterne:

\[ \hat{S}(t, \mathbf{X}) = \hat{S}(t)^{\exp(\mathbf{X}\hat{\beta})} \]

\(\hat{S}(t)\) er den “kontrollerede” \(S(t)\), der bestemmer overlevelsesforløbet for en person med specifikke karakteristika, hvor effekten af disse karakteristika er “kontrolleret” for.

Cox i praksis: baseret på flytte.csv.

library(tidyverse)
library(survival)
library(ggsurvfit)

df <- read_csv(here::here("flytte.csv"))

str(df)
spc_tbl_ [553 × 5] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
 $ nr     : num [1:553] 1 2 3 4 5 6 7 8 9 10 ...
 $ dur    : num [1:553] 12 16 19 22 26 27 28 29 32 33 ...
 $ status : num [1:553] 1 1 1 1 1 1 1 1 1 1 ...
 $ aargang: num [1:553] 6 6 6 6 6 6 6 6 6 6 ...
 $ udd    : num [1:553] 0 0 0 0 0 0 0 0 0 0 ...
 - attr(*, "spec")=
  .. cols(
  ..   nr = col_double(),
  ..   dur = col_double(),
  ..   status = col_double(),
  ..   aargang = col_double(),
  ..   udd = col_double()
  .. )
 - attr(*, "problems")=<externalptr>

Cox i praksis

library(tidyverse)
library(survival)
library(ggsurvfit)

df <- read_csv(here::here("flytte.csv"))

cox.m1 <- coxph(Surv(dur, status) ~ udd + factor(aargang), method = "efron", data = df)

cox.m1
Call:
coxph(formula = Surv(dur, status) ~ udd + factor(aargang), data = df, 
    method = "efron")

                     coef exp(coef) se(coef)      z        p
udd               -0.7314    0.4813   0.2195 -3.332 0.000863
factor(aargang)7  -0.1493    0.8613   0.2654 -0.562 0.573789
factor(aargang)8  -0.5364    0.5848   0.2490 -2.154 0.031229
factor(aargang)9  -0.1917    0.8255   0.2394 -0.801 0.423208
factor(aargang)10 -0.5713    0.5648   0.2706 -2.111 0.034785
factor(aargang)11 -0.7677    0.4641   0.3271 -2.347 0.018910
factor(aargang)12 -0.9039    0.4050   0.3368 -2.684 0.007277
factor(aargang)13 -1.0231    0.3595   0.3054 -3.350 0.000807
factor(aargang)14 -0.7809    0.4580   0.2902 -2.691 0.007133
factor(aargang)15 -0.5052    0.6034   0.3049 -1.657 0.097477
factor(aargang)16 -0.2959    0.7439   0.2830 -1.045 0.295883
factor(aargang)17 -0.9296    0.3947   0.2852 -3.260 0.001116
factor(aargang)18 -0.7148    0.4893   0.2607 -2.741 0.006122
factor(aargang)19 -0.4224    0.6555   0.3230 -1.308 0.190993
factor(aargang)20 -0.6113    0.5426   0.2815 -2.172 0.029888
factor(aargang)21 -0.5190    0.5951   0.5049 -1.028 0.304071
factor(aargang)22 -0.5530    0.5752   0.5371 -1.030 0.303241

Likelihood ratio test=63.42  on 17 df, p=2.835e-07
n= 553, number of events= 434

Hvad betyder proportionalitet ift. dette output?

  • Vi kommer senere til at se nærmere på hvordan vi tester antagelsen om proportionalitet.

Cox i praksis: faktor variable og ratios

Hvorfor vil vi gerne have kategorien med den største risiko som reference?

02:00

library(tidyverse)
library(survival)
library(ggsurvfit)

df <- read_csv(here::here("flytte.csv"))

cox.m1 <- coxph(Surv(dur, status) ~ fct_rev(factor(udd)) + fct_rev(factor(aargang)), method = "efron", data = df)

cox.m1
Call:
coxph(formula = Surv(dur, status) ~ fct_rev(factor(udd)) + fct_rev(factor(aargang)), 
    data = df, method = "efron")

                               coef exp(coef) se(coef)      z        p
fct_rev(factor(udd))0       0.73136   2.07790  0.21952  3.332 0.000863
fct_rev(factor(aargang))21  0.03404   1.03462  0.60585  0.056 0.955200
fct_rev(factor(aargang))20 -0.05831   0.94336  0.51210 -0.114 0.909351
fct_rev(factor(aargang))19  0.13063   1.13955  0.50602  0.258 0.796289
fct_rev(factor(aargang))18 -0.16177   0.85063  0.51066 -0.317 0.751402
fct_rev(factor(aargang))17 -0.37666   0.68615  0.51212 -0.735 0.462046
fct_rev(factor(aargang))16  0.25712   1.29321  0.53721  0.479 0.632202
fct_rev(factor(aargang))15  0.04774   1.04890  0.54839  0.087 0.930623
fct_rev(factor(aargang))14 -0.22788   0.79622  0.54020 -0.422 0.673144
fct_rev(factor(aargang))13 -0.47013   0.62492  0.54842 -0.857 0.391311
fct_rev(factor(aargang))12 -0.35089   0.70406  0.56710 -0.619 0.536083
fct_rev(factor(aargang))11 -0.21473   0.80676  0.56124 -0.383 0.702019
fct_rev(factor(aargang))10 -0.01831   0.98186  0.53059 -0.035 0.972471
fct_rev(factor(aargang))9   0.36127   1.43515  0.51561  0.701 0.483507
fct_rev(factor(aargang))8   0.01658   1.01671  0.51964  0.032 0.974552
fct_rev(factor(aargang))7   0.40372   1.49738  0.52830  0.764 0.444759
fct_rev(factor(aargang))6   0.55299   1.73844  0.53714  1.030 0.303241

Likelihood ratio test=63.42  on 17 df, p=2.835e-07
n= 553, number of events= 434

Cox i praksis

library(tidyverse)
library(survival)
library(ggsurvfit)

df <- read_csv(here::here("flytte.csv"))

cox.m1 <- coxph(Surv(dur, status) ~ udd + aargang, method = "efron", data = df)

cox.m1 %>% 
  survfit2() %>%
  ggsurvfit() +
  add_confidence_interval()

Cox i praksis

library(tidyverse)
library(survival)
library(ggsurvfit)

df <- read_csv(here::here("flytte.csv"))

km.m1 <- survfit(Surv(dur, status) ~ 1, type = "kaplan-meier", data = df)

km.m1 %>% 
  ggsurvfit(
    color = "red") +
  add_confidence_interval()

Cox i praksis (sammenlignet med KM)

library(tidyverse)
library(survival)
library(ggsurvfit)

df <- read_csv(here::here("flytte.csv"))

km.m1 <- survfit(Surv(dur, status) ~ udd, type = "kaplan-meier", data = df)

km.m1 %>% 
  ggsurvfit() +
  add_confidence_interval()

Cox i praksis (sammenlignet med KM)

library(tidyverse)
library(survival)
library(ggsurvfit)

df <- read_csv(here::here("flytte.csv"))

km.m1 <- survfit(Surv(dur, status) ~ udd + aargang, type = "kaplan-meier", data = df)

km.m1 %>% 
  ggsurvfit() +
  add_confidence_interval()

Cox i praksis

Hvis vi vil vide hvordan overlevelse varierer i henhold til en specifik covariat, kan vi stratificere analysen og stadig “kontrollere” for baggrundsvariable. \(S(t)\) estimeres for hvert strata, som her er uddannelsesretning, kontrolleret for årgang.

library(tidyverse)
library(survival)
library(ggsurvfit)

df <- read_csv(here::here("flytte.csv"))

cox.m2 <- coxph(Surv(dur, status) ~ strata(udd) + aargang, method = "efron", data = df)

cox.m2 %>% 
  survfit2() %>%
  ggsurvfit() +
  add_confidence_interval()

*_join() som forudsætning for dagens øvelse.

Hvor bekendte er i med at “joine” data? Skal vi køre en introduktion inden øvelserne?

Dagens øvelse

Vi kigger nærmere på alder for fødsel af første barn (inklusiv alderen på faren ved fødsel).

load(".../forlob_06032023.rda")

Datasættet minder om øvelsesdata fra sidste gang.

Det adskiller sig ved at det

  1. er person-period data (indeholder alle obs. frem til hændelse eller datastop)

  2. hændelsesvariable skal dannes ud fra C_ANTBOERNF

  3. regionsvariablen er nu tidsvaerierende

Datamanagement

  1. Lav region (tidsvarierende variabel) om til en tidsinvariant variabel, der indikerer regionsbopælen i personens 18. leveår (forløbets start). Den skal ligne region18 fra sidste lektion. hint: *_join()

  2. Konstruer ud fra C_ANTBOERNF en event-variabel: Tag højde for hvordan C_ANTBOERNF er defineret af DST (læs dokumentation på DST Times… Google)

    • Diskuter hvordan defininationen kan være et problem og hvordan vi kan løse problemt… hvis det er et problem?

  3. Dan en variabel, der er en indentikator for den sidste observation i forløbet. Den kan enten være en hændelse eller en censurering.

  4. Reducér data til kun at indeholde den sidste observation i forløbet.

Analyse

  1. Lav en KM
    • Forskellige tests - analysér og diskutér
    • Lav et KM plot
  2. Cox-regression (bivariat)
    • Bivariat med køn - analysér og diskutér hvad resultatet betyder
  3. Cox-regression (multivariat)
    • Indsæt køn, region (tidsinvariant) og uddannelse
    • Fortolk resultaterne
    • Plot kønnenes overlevelse med udgangspunkt i cox-regressionen

Jeppe F. Larsen | 6. marts 2025